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We consider a system of two coupled Tomonaga-Luttinger liquids (TLL) on parallel chains and 
study the Renyi entanglement entropy S„ between the two chains. Here the entanglement cut is 
introduced between the chains, not along the perpendicular direction as used in previous studies 
of one-dimensional systems. The limit n — ► 1 corresponds to the von Neumann entanglement 
entropy. The system is effectively described by two-component bosonic field theory with different 
TLL parameters in the symmetric/antisymmetric channels as far as the coupled system remains in 
a gapless phase. We argue that in this system, S„ is a linear function of the length of the chains 
(boundary law) followed by a universal subleading constant 7„ determined by the ratio of the two 
TLL parameters. The formulae of y n for integer n > 2 are derived using (a) ground-state wave 
functionals of TLLs and (b) boundary conformal field theory, which lead to the same result. These 
predictions are checked in a numerical diagonalization analysis of a hard-core bosonic model on a 
ladder. Although the analytic continuation of y n to n — > 1 turns out to be a difficult problem, our 
numerical result suggests that the subleading constant in the von Neumann entropy is also universal. 
Our results may provide useful characterization of inherently anisotropic quantum phases such as 
the sliding Luttinger liquid phase via qualitatively different behaviors of the entanglement entropy 
with the entanglement partitions along different directions. 

PACS numbers: 71.10.Pm, 03.67.Mn, 11.25.Hf 



I. INTRODUCTION 

The concept of Tomonaga-Luttinger liquid (TLL) pro- 
vides a universal framework for studying various one- 
dimensional (ID) interacting systems .El ThoJow- lying ex* 
citations of such a system, either fermioniccl or bosonic,H 
are essentially collective, and can be recast into a bosonic 
field theory describing the density and phase fluctua- 
tions. A spinless TLL is characterized by a continuously 
varying parameter K (so-called TLL parameter), which 
appears in the exponents of correlation functions in the 
ground state and experimentally in the power-law tem- 
perature dependence of response functions. When two 
spinless TLLs are coupled (or when an interaction is in- 
troduced in a ID gas of spin-i particles), the bosonic 
fields are reorganized into symmetric and antisymmetric 
channels, which can independently form TLLs. This is 
a fundamental mechanism which also underlies the spin- 
charge separation in a ID electron gas. Interestingly, 
this idea has been generalized to a two-dimensional (2D) 
array of coupled TLLs, predicting a novel non-Fermi liq- 
uid phase, called sliding Luttingcr liquid, which shows 
highly anisotropic correlations.ui3 A fundamental ques- 
tion related to these studies is in what way the system of 
coupled TLLs are distinguished from more conventional 
phases such as Fermi liquids or from the decoupled TLLs. 
Stimulated by the recent advances in applying quantum 
information tools to many-body systems, we here address 
this question using one of such tools — the entanglement 
entropy in the ground-state wave function. 

By partitioning the system into a subregion A and its 
complement A, the entanglement entropy is defined as 
the von Neumann entropy Sa = — Trp^ In pa of the re- 



duced density matrix pa = Tr^|4')(^ r |, where |\&) is the 
ground state of the system. When the system contains 
only short-range correlations, A and A correlate only 
in the vicinity of the boundary, and the entanglement 
entropy Wfflcs with the size of the boundary (bound- 
ary law) Del Deviation from the boundary law signals the 
presence of certain non-trivial correlations, and further- 
more can contain universal numbers characterizing the 
system. In one-dimensional critical systems, for example, 
the entanglement entropy Sa for an interval embedded in 
the system shows a logarithmic scaling, whose coefficient 
reveals the central Cjkacge c of the underlying conformal 
field theory (CFT).etl3 Possible further information of 
CFT such as the TLL parameter K caia-be encoded in a 
multi-interval entanglement entr.ppvEjlLj and in correc- 
tions to the universal scalings.EjEH. In topologicall^-ns* 
dered systems&E3 and in some 2D critical systems ,EIrEI3 
the entanglement entropy obeys a boundary law, but 
there appears a subleading universal constant which is 
determined from the basic properties of the ground state. 

In this paper, we aim to characterize the quantum en- 
tanglement arising from the coupling of TLLs. We con- 
sider, as the simplest situation, a system of two coupled 
spinless TLLs defined on parallel periodic chains (rings), 
and study the entanglement entropy between the two 
chains. The system is effectively described by a two- 
component bosonic field theory with different TLL pa- 
rameters K± in the symmetric and antisymmetric chan- 
nels. If we identify the two chains with the spin-i degrees 
of freedom, these channels correspond to the charge and 
spin modes, respectively. For ID systems, the entangle- 
ment entropy has so far been studied mostly for an inter- 
val embedded in the chain, which can count the central 
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charge in critical systems. We here instead partition the 
system into two rings. We expect that this partitioning 
is more useful in observing the effects of the coupling of 
the two TLLs. Furthermore, we expect that the present 
setting provides a good starting point for understanding 
possibly highly anisotropic characters of entanglement in 
a 2D sliding Luttingcr liquid. 

Specifically, we construct the reduced density matrix 
Pa for one of the chains by tracing out the other, and 
compute the Renyi entanglement entropy: 
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n - 1 



log(Tr p\) 



(1) 



The limit n — > 1 corresponds to the von Neumann entan- 
glement entropy: 



Si = lim S n = -Tr p A \ogp A - 



(2) 



The limit n — j'-yo corresponds to the so-called single-copy 
entanglemcntO: 



log Ar, 



(3) 



where A ma x is the largest eigenvalue of pa- It is known 
that the values of S n with integer n > 2 determine 
the full eigenvalue-distribution of pa (so-called entangle- 
ment spectrum) £3 When there is no coupling between 
the chains, S n is simply equal to zero. The entropy S n 
increases as the coupling increases. 

We will see that S n (with 11 = 1,2,..., 00) obeys a 
linear function of the chain length L: 



S n = a n L + 7„ 



(4) 



where the ellipsis represents terms which are negligible 
in the limit L — > 00. The first term a n L can be simply 
viewed as a boundary law contribution, and the coeffi- 
cient a n depends on microscopic details. Our main in- 
terest lies in the subleading constant 7„. We argue that 
this constant is universal and is determined by the ratio 
of two TLL parameters- K + / A'_ . 

Recently, PoilblancEa studied the entanglement en- 
tropy for a similar partitioning in gapped phases of a spin 
ladder model. In his results, the entanglement entropy 
shows a similar linear scaling, but a subleading constant 
was not identified. We expect that the linear scaling is a 
generic feature of this type of partitioning, and that the 
appearance of the subleading constant is characteristic of 
critical systems. 

The paper is organized as follows. In Sec. |n|, we set up 
the problem which we consider in this paper, and present 
path integral representations of the reduced density ma- 
trix moments Tr p n A (with integer n > 2). Based on these 



representations, in Sees. [II and IV, we calculate the mo 



ments using two different approaches. In Sec. [II, we use 
the field theoretical representations of the TLL ground- 
state wave functions. In Sec. IV, we use a modern tech- 



compactification lattices. The two approaches are com- 
plementary: while the derivation is simpler in the former, 
the latter provides a more systematic treatment which 
does not require any regularization procedure. Both the 
approaches lead to the linear scaling of S n and the same 
formulae for the subleading constant 7„. The expres- 
sions of 7„ ( as a function of K+/K-) are summarized 



in Sec. [II C, together with a discussion on their ana- 
lytic properties. In particular, we discuss a difficulty in 
analytically continuing the formulae of 7„ (obtained for 
integer n > 2) to the von Neumann limit n — > 1 + . In 
Sec. 0, we check our predictions on S n (n > 2) in a nu- 
merical diagonalization analysis of a hard-core bosonic 
model on a ladder. While we do not have any analytic 
prediction on the von Neumann entropy Si, the numer- 
ical result suggests that Si obeys a linear scaling as the 
Renyi entropies does and that the subleading constant 71 
is universal. We conclude with a summary in Sec. VI. Im- 



plications of our results on a 2D sliding Luttinger liquid 
are also presented. 



II. SETUP OF THE PROBLEM 

In this section, we set up the system and the problem 
which we consider in this paper. In particular, we present 
the path integral representations of the reduced density 
matrix moments Tr p\ (with integer n > 2), which will 
be used in the following sections. 



A. Coupled Tomonaga-Luttinger liquids 

We consider a system of two TLLs H u (v = 1, 2) on 
parallel periodic chains of length L coupled via interac- 
tions Hi2- We assume that the two TLLs are equivalent 
and are described by the Gaussian Hamiltonian: 
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1,2, 



nique in boundary CFT based on boundary states and 



(5) 

where x is the coordinate along the chains, and v and 
K are the velocity and the TLL parameter, respectively, 
in each chain. In the case of fermions, K < 1 (K > 1) 
corresponds to a repulsive (attractive) intra-chain inter- 
action. The dual pair of bosonic fields, <fi v and Q v , sat- 
isfy [d?„{x),9 v >{x')] = (i/2)[l + sgn(a; - x')]5 vvl . The 
field 4>v( x ) is related to the particle density p u {x) via 
Pv(x) Ri po — -^7^^%^ with po being the density in the 
ground state while the field Q v (x) represents the Joseph- 
son phase. We assume that there is no particle tunneling 
between the chains, and therefore the particle number is 
separately conserved in each chain (f/(l) x £7(1) symme- 
try). 

Now let us consider, for instance, the interaction of the 
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form 



H 



12 



U d(j>i d<j>2 



i 7r dx dx 



(6) 



which corresponds to the leading part in the density- 
density interaction. To treat this, we introduce the sym- 
metric/antisymmetric combinations of the bosonic fields: 



±02), 



>i ± fe). (7) 



Then the total Hamiltonian H = Hi + H-2 + -H12 can be 
formally decoupled into two free bosons defined for these 
symmetric/antisymmetric channels: 



with 



H . 



v± 



dx 
2 



H = H-l 



A' 4 



H- 













\ dx ) 



(8) 



(9) 



Here the renormalised velocities v± and TLL parameters 
K± are given by 



v± 



1 ± 



KU 



K+ = K 1 ± 



A[A 

7TW / 



(10) 



Note that, although the two channels are formally de- 
coupled in Eq. (@), zero modes of the two channels are 
intertwined, which will be seriously discussed in Sec. IV. 
On the other hand, the oscillator modes of these channels 
are completely decoupled. In general, if if 12 consists only 
of forward scattering processes, the total Hamiltonian if 
can be similarly recast into the form in Eqs. (||) and (||). 
Even when H12 contains other terms, this form is still 
applicable as long as those terms are irrelevant and di- 
minish to zero in the rcnormalization group (RG) flow. 
In this case, K± can change slightly from the perturba- 
tive result [like Eq. (|l0|)1 along the RG flow, and their 
precise values in the infra-red limit can be determined 
by examining correlation functions numerically, for ex- 
ample. In the following, we consider the situation where 
the Hamiltonian in Eqs. (||) and @ presents the exact 
long-distance physics, and treat K± as free parameters. 
Since we are interested in the entanglement between the 
two chains, we keep in mind that the bosonic fields <j)± 
and 9± diagonalizing if are related to the original fields 
on the chains via Eq. (^). Note that Eq. (0) is protected 
by the permutation symmetry of the two chains and is 
applicable beyond the perturbative regime. 



B. Path integral representations of reduced density 
matrix moments 

For the ground state |^) of if, we consider the Rcnyi 
entanglement entropy S n [Eq. (Q)] with integer n > 2 
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FIG. 1: The path integral representation for the reduced 
density matrix pA in Eq. (|l5|). In each of the two sheets, 
a periodic boundary condition is imposed in the x direction. 
Therefore the left and right sheets form a cylinder and a torus, 
respectively. 



between the two chains. Here we represent the moments 
of the reduced density matrix, Tr p n A , in the language of 
the path integral. We start from the finite-temperature 
density matrix of the total system: 



1 



-/3H 



with Z = Tr 



-8H 



(11) 



The inverse temperature (3 is eventually taken to infinity 
so that p —} \^)(^f\. We move on to the path integral 
formalism in the Euclidean space time (tE,x). The Eu- 
clidean action is 



Se = 



dt, 



[ dx (£ 
Jo 



with 



C 



E± 



e++C E -) (12) 



4±) 2 +v± 2 (dtM 2 ] ■ (13) 



Although Se is diagonalizcd in <fi± basis, in the follow- 
ing, we rather regard this as a functional of </> 12 using 
the relation (f?J) since we are interested in the entangle- 
ment between the two chains. On this ground, the matrix 
element of the density matrix p is expressed as 



1 



V„(*) 



T>(f>iT>4>2 e 



— Se [<t>i,<t>2] 



(14) 



where ip v = {^jy(a;)}o<2;<L and those with a prime are 
field configurations defined along the chains 1 and 2 re- 
spectively. The path integral is done under the condition 
that 4> v {i,E,x) (y = 1,2) is equal to <f v {x) and <p' u (x) at 
the imaginary time = and /3, respectively. 

The reduced density matrix pa for the chain 1 is ob- 
tained by identifying if 2 and ip' 2 in Eq. ( fTi] ) and integrat- 
ing over (p2~. 



{<Pi\Pa\<Pi) : 

<1>i(J3,x) 



Vip 2 (ip[,ip2\p\(pi,(p2) 



= <Pi(x) vfavfa e 

= f'i(x) 



— Se[01,02] 



(15) 
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FIG. 2: The path integral representation of Trp^ in Eqs. jlrj ) 
and (|l^), for the case of n — 3. 



with T>tp v = Y[ x d<£ u {x). We introduce a graphical repre- 
sentation in Fig. [j], where two sheets express the space- 
time on which the fields 0x,2 are defined. The partial 
trace in Eq. ( |l5| ) corresponds to gluing the two edges of 
the sheet for (/> 2 . 

Now we consider the n-th moment of the reduced den- 
sity matrix, Trp^, with integer n > 2. To construct this, 
we consider n copies of the diagram in Fig. |l| and glue 
them cyclically, as illustrated in Fig || for the case of 
n = 3. This leads to an expression 



/n n 
JJXVaj-l Y[{<P2j+l\pA\<P2j-l) 
3=1 3=1 

/2n n 
n p ^ n ( < ^ 2 J+ i ' ^ i/'i^-i ' 
„- 1 „' 1 



(16) 



3=1 3=1 



where 's with odd (even) subscripts are defined for the 
chain 1 (2) and <^2n+i = fi- This can be expressed in a 
compact way: 



~z<> 



(17) 



where Z n is the partition function defined for 2n sheets 
which are interconnected as shown in Fig. |^. The dia- 
gram consists of a large torus for </>i's, and n small tori 
for fa's. Because of the interactions between and 02, 
the calculation of such a partition function is not trivial. 
In Sees. [II and |lV|, we present different ways to compute 
Eq. ([if]) or Eq. (p-7[), which eventually lead to the same 
result. Here we mention the case of no inter-chain inter- 
action H\2 = 0, where the tori in Fig. || are decoupled. 
Using the ground-state energy Eq of H\ and H2 , we have 
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FIG. 3: Rewriting of Z3 in Fig. using cj>± basis. The par- 
tition function Z3 is obtained after integrating over the field 
configurations tpi, . . . ,(p6- 



Z n « e~ n ? E <> (e-^ )" and Z n « (e" 2 ^)" in the limit 
j3 — > 00, which lead to S„ — 0. 



III. WAVE FUNCTIONAL APPROACH 

In this section we compute Eq. ( |l6| ) using a field the- 
oretical representation of the TLL wave function, and 
derive the expressions of the Rcnyi entropies S n for inte- 
ger n > 2. Similar approaches were also used to calculate 
the entanglement entropy in 2D criticaljjraizje functionsEj 
and the ground-state fidelity in TLLs.l2tl3 In this sec- 
tion, we do not include the zero modes of the bosonic 
fields and regard H± as completely independent. This is 
justified because we arc interested in the entanglement 
properties of the ground state, where zero modes do not 
appear. 



Reduced density matrix moments and wave 
functionals 



The difficulty in computing Z n comes from the inter- 
actions between different sheets in Fig. ||. To treat these 
interactions, we work in the symmetric/antisymmetric 
basis of bosonic fields, in which the action is diagonal, 
leading to a diagram as in Fig. | As a tradeoff, the 
boundaries of the sheets are now interconnected in a non- 
trivial way. The strategy of this section is to first treat 
each sheet of Fig. |3| separately by fixing the boundary 
field configurations, (f/s, and to then integrate over y?,-'s 
to calculate the partition function Z n . 
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As mentioned in Sec. [II A , the winding numbers (zero 
modes) of the symmetric/antisymmetric channels are in- 
tertwined, and therefore the two channels are not com- 
pletely decoupled. However, since we are interested in the 
entanglement properties of the ground state in the limit 
(3 — > oo, we can work in the sector of the Hilbert space 
where the winding numbers are set to zero. Namely, 
we focus on the oscillator modes in the Hamiltonian. 
In this sector, H± commute with each other. Using 
e~@ H = e~ l3H +e~ l3H - , we rewrite the matrix element of 
p appearing in Eq. (|T|) as 



(V>2j+l,V2j|/#2i-l,V>2i) 
= 1_ I ip 2 j+l + P2j p -fj H ^ 

~ Z\ V2 

I <f2j + l - f2j 



A V2 



-PH- 



<P2j-l + <P2j \ 

V2 I (18) 

(p 2 j-l - <fi 2 j \ 

V2 r 



An expression of the form ((p'\e~^ H± \ip) in this equation 
corresponds to each sheet in Fig. |^. Since H± are the 
Hamiltonians of massless free bosons, (p'\e~^ H± \p) can 
be viewed as the propagator of a closed bosonic string in 
the imaginary time. Such a "closed string propagator" 
has been computed in, e.g., Refs. [33] [55] [in particular 
a compact expression is shown in Eq. (24) of Ref. [55| . 
Rather than using the expression obtained in these works, 
we here take a simpler route as follows. 

We take the limit p — > oo, and then only the ground 
states |$±) of H± (with eigencnergies E±) contribute to 
the propagators and the partition function: 



Z = Tier 13 " ss e -P(E++E-)_ 



(19) 
(20) 



Using these, we can rewrite Eq. ( |T8| ) as 

(ip2j+i,(p2j\p\<P2j-i,V2j) 



V2j 



lf2j+l - <P2j 



V2 

¥>2j-l - V2j \ 



/ (21) 



V2 /■ 



Here, an expression of the form (p\^±) is the represen- 
tation of a ground state wave function in terms of the 
field configuration {<p(x)}o< x <L along the chain, which 
we call a "wave functional" following Ref. Bq. 



B. Calculation of reduced density matrix moments 

The ground state |-wpjp-fui3ctional of a TLL has been 
derived in literature .E3cilc3n2j In Appendix |A], we present 
its si mple derivation in the operator formalism. From 
Eq. flAl^ ), the wave functional is expressed as 



£[¥>] 



(22) 



where £[ip] is a quadratic functi onal of ip [see Eq. (A.14) 
for the explicit form]. From Eq. (All), the normalization 
factors Af± are given by 



A4=n 



(23) 



Using Eqs. © and © and £[ip] = £[-ip], Eq. © is 
rewritten as 



. 2n 



x exp 



2n , 

El 



ipj +tfj+i 



V2 



K- 



V2 



(24) 

Here each term in the argument of the exponential func- 
tion corresponds to an edge of a sheet in Fig. ||, and 
represents the probability distribution of field configura- 
tions in a wa y an alogous to the Boltzmann weight. As 
shown in Eq. (A9), the functional £[p] has a very simple 
form when written in terms of the Fourier components 
{ip m } of p: 



£[{0m}} = k 



m | *Pm 



(25) 



Therefore, we ex pan d pj into the Fourier components 
{<Pj, m } as in Eq. QA6]), and rewrite Eq. (24) as 

„ 2n oo 

Tr P \ =(N+N-r n / n n (<%™#i,m) 

** j=lm=l 
/ oo 

x exp I - — ,„ $t,M 



where <& m = (tpi, m ,tp2, m , ■ ■ ■ , (p2n > m) t and M n is a 2nx2n 
matrix defined as 




(26) 



(A \B 
\B A \B 



\B A 



(27) 



\B A J 




Performing the Gaussian integration and using Eq. (|23j) , 
Eq. ( p6| ) is calculated as 

Tr p'X = (M+M-)- n ] 



m— 1 



l 



detM„ 



[ (dctM„)' 1 



(29) 



G 



Since M n has the same form as the Hamiltonian of a ID 
tight-binding model, it can be easily diagonalized and its 
determinant is calculated as 



2n-l 



1=0 



det M n = TT A;, A; := A + Bcos ( ^ 

\ 2n 



(30) 



In Eq. we have obtained an infinite product of 

the form []^ =1 C _1 (with C > 1), which needs to be 
regularized. We introduce a short-distance cutoff clq of 
the order of the lattice spacing. We rewrite the prod- 
uct as Wjn^C" 1 / 2 . In this expression, m runs over 
L/ao — 1 modes by considering the exclusion of the zero 
mode. Therefore the product scale as C 1 / 2 e _Qi (with 
a = (logC)/(2a ) > 0). The prcfactor C 1 / 2 gives a 
cutoff-independent (and thus universal) constant. [A 
similar technique has also been used for evaluating the 
fidelity in a TLL in Ref. [32|] We note that the same uni- 
versal constant can also be obtained by the ^-function 
regularization. Applying this argument to Eq. (p9|), we 
arrive at 



-aL 



(dot AT, 



\l/2 



(31) 



where a is a cutoff-dependent constant. Note that, al- 
though we initially assumed integer n > 2, the final 
expression ([}]]) also contains the case of n = 1, where 
Tr pa — 1. This can be seen by setting a — and 



M 2 



A B 
B A 



(32) 



Compared to Eq. (^7|), we have B instead of \B in the 
elements because the elements on the subdiagonal parts 
and at the upper-right/lower-left corners in Eq. (^7|) are 
combined. The determinant and the eigenvalues of M% 
are written in the same ways as Eq. (|30"|). 



C. Expressions of Renyi entropies 

Equation (^lj) leads to a linear scaling of S n as a func- 
tion of the chain length L as in Eq. (|]) . The coefficient a n 
of the linear term depends on the short-distance cutoff ao 
and therefore is not universal. The sublcading constant 
term j n for integer n > 2 is obtained as 



-1 



log(detM„) = 1 J2 lo § A <- ( 33 ) 
2(71 — 1 * — ' 



2n-l 



In = 



2(n-l) 



2(n - 1) 



We see that 7„ is determined by the underlying field the- 
ory and is a function of the ratio of the two TLL param- 
eters, K+/K-. As an example, for n — 2, one obtains 



72 = - log 




(34) 



In the limit of n — > oo, the summation over Z in Eq. (|3 
is replaced by an integral, leading to 



7oo 



with 



I(s) 



log 




(35) 



cZ6» 
2^ 



log(l 



log 



i + vr 



(36) 

Here the integral was calculated as follows. We differen- 
tiate I(s) with respect to s and integrate over 9: 



dl(s) 
ds 



2jv 



2tt 1 + s cos 9 



1 



3 VI 



(37) 



Using 1(0) — and integrating this over the interval [0, s] 
give the final expression in Eq. (|3^). 

In the replica procedure for calculating the von Neu- 
mann entropy, we compute the Renyi entropies S n for 
integer n > 2, take the analytic continuation to real 
n € [1, oo], and then take the limit n — > 1 + . In Eq. (|3^), 
we cannot find any obvious way to extend the formula 
of 7„ to the case of real n. Let us focus on the expres- 
sion 7„ := log(detM n ) appearing in Eq. (|33|). To gain 
insights on the analyticity of j n as a function of n, we 
expand it around K + /K^ = 1, which corresponds to the 
limit of no inter-chain coupling. To this end, it is useful 
to introduce a parameter 



- K. 



K- + If- 



Using this, 7„ is written as 



dog(l- K 2 )+ Y, log 



1=0 



1 + K COS 



2tW 
~2n 



Expanding around k — gives 



A nm 2m 
K 

m—1 



2m 



(38) 



(39) 



(40) 



with 



1=0 v 7 

1 2771 / v In— 1 

1 yf 2 "M v p 2 " 

fe=o v 7 ;=o 



(41) 



l(k—rn) /n 



In the summation over k, only the terms where fc — to is 
an integer multiple of n contribute. For m < n, it occurs 
only for k = to, and A n/m is given by a simple expression 



2n /2m 
^»,m - 22m I m 



(42) 
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For m > n, A n m can contain other terms and show non- 



trivial dependences on n and m. For example, for m 
one obtains 

A -I 2 {n = 1) 
" ,1_ ln (n>2). 



1, 



(43) 



This leads to the lowest-order expansion of 7 n for n > 2 



5« 2 + o(« 4 ), 



(44) 



which is not smoothly connected to 71 = as n — > 1 
Multiplying -l/[2(n - 1)] to Eq. (El), we obtain 



4(n- 1) 



(45) 



for n > 2. If we naively take the limit n — > 1 + in this ex- 
pression, we find that the coefficient of the leading (ordcr- 
k 2 ) term in j n is divergent. This indicates that in this 
problem, it is not easy to study the von Neumann en- 
tropy Si from the knowledge of the Renyi entropies S n 
with integer n > 2. At present, we do not have any ana- 
lytic prediction on Si. However, our numerical result in 
Sec. [v] indicates that Si also obeys a linear function of 
L and that the subleading constant 71 is determined by 
K+/K-. In particular, for small k, wc find that 71 obeys 
a non-trivial power function 



7i 



(46) 



with 6 ps 1.6-1.7. In spite of the qualitative difference 
between Eq. ( fi"5|) and Eq. (E^) , our numerical result also 
suggests that for fixed k, j n changes rather smoothly 
when n is changed from 2 to 1. In Sec. [y], wc will present 
some possible scenarios as to how these two different 
small-K behaviors arc connected to each other. 

The current problem adds to the list of problems where 
the Renyi entropy shows quite a non-trivial analyticity as 
a function of n. Here we cite a few examples known in. 
literature. In massive integrable quantum field theory,c2l 
the analytic continuation could not be uniquely intro- 
duced from the knowledge for integer n > 2, and the 
appropriate one needed to be chosen carefully. In the 
Renyi-entanglement entropy of two disjoint intervals in 
CFT,t3 the analytic form for integer n > 2 has a non- 
trivial form, and its analytic continuation to n — > 1 + has 
been achieved in certain limits, leaving a general solution 
open. In thqJlcnyi entropy of a line embedded in 2D 
Ising models,c3 the constant part behaves as a step-like 
function of n with discontinuity at n = 1, which means 
that the standard replica procedure fails to address the 
case of n = 1. 



IV. BOUNDARY CONFORMAL FIELD 
THEORY APPROACH 

In this section we express the partition functions, Z n 
and Z n , in Eq. ( p"7j ) as the transition amplitudes between 



conformal boundary states. In the limit j3 L ^> 1, 
these partition functions contain universal multiplicative 
constant contributiops, known as the boundary "ground- 
state degeneracies. "cil This approach does not require 
any regularization procedure and determines the univer- 
sal contributions in the partition functions in a way con- 
sistent with a certain condition under the_modular trans- 
formation (Cardy's consistency conditioner) . Similar ap- 
proaches were also used quite recently to calculate—!, 
entanglement entropy in 2D critical wave functionscj 
and the ground-state fidelity in TLLs.e3 



A. Compact ificat ion conditions of bosonic fields 

To app ly boundary CFT to the system introduced in 
Sec. II A, one needs to precisely discuss the compactifica- 
tion conditions imposed on the bosonic fields. Wc will sec 
that zero modes of the symmetric/antisymmetric chan- 
nels are intertwined and require a careful treatment. 

The original bosonic fields, 4> v and (y = 1,2), de- 
fined along the chains are compactified on circles with dif- 
ferent radii. When periodic boundary conditions (PBC) 
are imposed, these fields can acquire aanding numbers 



when going around the chains, namely,cij 

<j) v {L) = <j) v (Q))+2irrn v , 
OAL) = MO) + 1-Kfm v , 
n„, m v G Z. 

Here the compactification radii arc given by 



(47) 



2^' 



(48) 



Before discussing the compactification conditions in 
the symmetric/antisymmetric channels, let us mention 
that H in Eq. (||) is not a Hamiltonian of a conformally 
invariant system because the two velocities v± in Eq. (^) 
are different in general. To apply boundary CFT later, 
we restore the conformal invariance by simply replacing 
v± —> 1. Although this replacement changes the spec- 
trum of the Hamiltonian, it does not change the eigen- 
states. The ground state and therefore its entanglement 
properties should remain unchanged. To further sim- 
plify the Hamiltonian, we rescale the bosonic fields as 

$± = 0±/v~ 
nian is 



and 0± = JK±Q±. The new Hamilto- 



H = - 
2 



dx 




(49) 



with 



$ = 



e = 



©+ 



(50) 
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FIG. 4: Two representations of Zi. Four sheets are interconnected in (a). We fold each sheet at is = /3/2 and superpose all 
the 8 pieces of sheets, leading to a 8-component field $ on a single sheet in (b). 



From Eq. (pT7|), the new bosonic fields are subject to 
the conditions 

$(L) = $(0) + 2nu, u = nidi+n 2 a 2 , (51) 

9(L) = 9(0) + 2nv, v = m 1 b 1 +m 2 b 2 , (52) 
n 1 ,n 2 , mi, m 2 € Z, 



where 



Ol,2 



V2 \±i/Vk. 




Note that <ii • bj = ^Sij . Let A be the lattice of u defined 
by Eq. (|5l|), and let A* be its reciprocal lattice. Then v 
defined by Eq. ( f52] ) lives on ^ A* . The lattices of u and 
w introduced in this way are called the compactification 
lattices of $ and 9. 



B. Reduced density matrix moments and 
boundary states 

We consider the partition functions, Z n and Z n , ap- 
pearing in Eq. (P~T|) . In the following discussions, we 
mainly focus on the case n = 2; the generalization to 
arbitrary integer n > 1 is straightforward and will be 
done in Sec. IVD. From the argument in Sec. 



[IB, Z- 



is expressed by four sheets interconnected as shown in 
Fig. [I](a). Invoking the folding technique of Refs. p6|fi4| , 
we fold each sheet at = (i/2 and superpose all the 
8 pieces of sheets. As a result, we have a 8-component 
bosonic field $ living on a cylinder of lengths L and /3/2 
in the spatial and temporal directions respectively; see 



Fig. |^(b). The Hamiltonian H for this "system" is writ- 
ten in the same form as in Eq. ([49]), but now $ and 9 
consist of 8 components each: 



(D <fi(l) ft(2) ifi(2) «(3) ffi (3) ffi (4) $ (4)^ f 



$ = ($v ' , , $V , , $v , $ W , $ 



(54) 



e = (e« , eL 1 ' , , eL 2) , e<?> , eL 3) , ef , e w )* . (55) 

Here, the components of $ are related to (j>J s in 
Fig. |(a) as 



U) 



± 



(56) 



and e2 -) are defined as their dual counterparts. The 8- 
component fields are subject to the conditions 



$(L) = $(0) +27TM, UE~ 

9(L) = 9(0) + 2ttu, tjgS 



A 



2tt 



-A" 



(57) 
(58) 



The primitive vectors of the lattice S are a„ [v 
1.2: j = 1,2,3,4), each of which is defined by insert- 
ing a v into (2j — 1)- and (2j)-th elements and zeros into 
the others. Similarly, the primitive vectors of the lattice 
H are bv (y = 1,2; j = 1,2,3,4), defined likewise from 
t- 

At the two boundaries at = P/2 and 0, the following 
boundary conditions are imposed respectively: 



r 2 : <f>? j) 

U = 







4 2i) 








0i 2J+1) , 












2 ; 4 ] - 







(i = 1,2), (59) 



(60) 
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We will express these conditions using boundary states, 
\Ti) and |r 2 ). The partition functions we wish to calcu- 
late are expressed as the transition amplitudes between 
these states: 



Z 2 
Z 2 



Z 



rir 2 



-Tiri 



(ri|e-^|r 2 ) 

{T x \e~^\Y x ) 



(61) 
(62) 



C. Boundary state formalism 



Before considering the two boundary states |IY 2 ) 
in more detail, we discuss the construction of bound- 
ary states in a more general setting. The bound- 
ary CFT for multicoraDoaent bosons has been devel- 
oped in string ttaocjCJEJ and applied to condensed 
matter problems .E3tiil In particular, a "mixed" Dirich- 
let/Neumann boundary condition, which we focus on 
here, has been discussed in Refs. ^7],||l]. Such "mixed" 
conditions have recently been applied to the calcula- 
tion of the. entanglement entropy in 2D critical wave 
functionsEj Here wc review basic knowledge on the 
boundary CFT for multicomponcnt bosons, and dis- 
cuss how to construct the boundary state for a "mixed" 
Dirichlet/ Neumann condition. For further details, we re- 
fer the reader to, e.g., Rcfs_ p6p9| , |5l| (especially Ref. |26| 
for the present applications), which contain useful sum- 
maries of boundary CFT for multicomponent bosons. 
The main result of this subsection is the formula of the 
boundary "ground-state degeneracy" in Eq. (^2|), which 
is used later to calculate universal (non-extensive) con- 
stant contributions in partition functions. 

We consider a c-componcnt free boson defined by the 
Hamiltonian H in Eq. (fl9j). The system is placed on a 
cylinder like Fig. |](b) and we impose certain conformally 
invariant boundary conditions at both ends. Since the 
PBC is imposed in the x direction, the bosonic fields 
have the following mode expansions: 

27T - 

$(t,x) =$o + —(ux + vt) (63a) 



L 



E 



i 



ik, 



(x-t) 



h.c 



27T - 

Q(t,x) = e + — (vx 

Li 



ut) 



(63b) 



E 



i 



i V47 



ikjn (x-\-t) 



(x-t) 



h.c 



with k m = 2irm/L. The spectra of u and v form the 
lattices 5 and S = j^S*, respectively. We have included 
the dependence on the real time t, which help to see 
that (a^) represents a left (right) moving mode. The 
elements of vectors, which we label by j = 1,2, . . . , c, 
obey the commutation relations 



L ■ a if 



r R -Rt 

l a m,ji a m',j' 



(64) 
(65) 



Using the expansions (|63|), the Hamiltonian H is diago- 
nalized as 



H 



2tt 



n(u 2 + v 2 ) + ^ 



l m) 



C 

12 
(66) 

where the last term comes from the zero-point motions 
of oscillators (Casimir effect). The ground state |$) of 

H is given by the condition a™' |*) = u\%) = v\$>) = 
We can decompose Eq. 
as 



3h into the chiral components 



e(t,x) = $ L (x + )-$ R (x-), 
x± = t ± x. 



(67) 



with 



®l/r(x±) = -($ ± 6 ) + j (±u + v 



E 



1 



t L/R e ^k m x ± , h ^ 



(68) 



We now introduce a conformally invariant boundary 
condition T at the time t = 0. Boundary conformal in- 
variance implies that the momentum density operator 
Tl — Tr vanishes at the boundary. Here, T^i R {t,x) = 

T L / R (x±) — 2ir(d±<&) 2 (with d± :— d x± ) are the chiral 
components of the energy-momentum tensor. The con- 
formal boundary state |T) therefore satisfies 



[T L {x)-T R {x)]\T)=Q. 



(69) 



[Here, T L / R {x) is defined by T L / R (t — 0,x). The same 
convention applies to 3>l/_r(x) and Jl/ R (x) below.] In a 
multicomponent boson, one can also introduce additional 
symmetry requirement of the form 



Jl(x) - KJ R (x) \T) = 



(70) 



which represents the conservation of currents in a general 
form (associated with a Hciscnberg algebra). Here 1Z is 
an orthogonal matrix, and 

•h/ R (t, x) = J L / R (x±) = ±d ± @(t, x) = d ± $ L/R (x ± ) 

(71) 

are the chiral components of the current operator. Since 
T L/R (x±) = 2n(Ju R (x±)) 2 , Eq. © implies Eq. ©. 
Therefore, Eq. ( J70| ) defines a subclass of conformal 
boundary states for multicomponent bos pna, .phich have 
many interesting physical applications p^I-P 1 ! On the 
other hand, conformal boundary states which satisfy only 
Eq. ( |69| ) and not Eq. (|7(i| ) are also known.c3 

We now focus on the subclass defined by Eq. ( |70| ) . The 
condition can be rewritten as 



d x ® l {x)+k$ r {x) |r> = 



(72) 
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This means that <f>L + IZQr is fixed at a constant vector 
along the boundary. In particular, setting TZ to the iden- 
tity matrix /, we have the Dirichlct boundary condition 
("D"), where 4> is fixed at a constant vector along the 
boundary. Setting TZ = —I leads to fixing 6, which then 
means the Neumann boundary condition ("N") dt& = 
(since d t <£> — d x Q). 

To obtain the explicit form of we decompose 
Eq. ( f72| ) into Fourier components using Eq. (|68|), lead- 
ing to 

+ |r)=0, (73a) 

(ai + 1Za%)\T) = (i&t + 1Za? n )\T) = 0. (73b) 

The solution of Eq. ( |73b| ) is given by the Ishibashi statcll 
\(u, v))) := exp ( - J2 3% . TZafA \(u, v)), (74) 

V m=l / 

where \(u,v)) is an oscillator vacuum characterized by 
the zero mode quantum numbers (or "winding numbers" ) 
u G H and v G 3. If (u, u) satisfies 



(u + v) + K(u -v)=0 



(75) 



required from Eq. (73a), the Ishibashi state \ (u, v))) satis- 
fies the conformal invariancc. It is known, however, that 
in order to obtain a stable boundary state for a given 
TZ, one must take a linear combination of the Lshibashi 
states over all possible (u,v) satisfying Eq. (|75|)£j 

We proceed our discussion focusing on the case of a 
"mixed" Dirichlet /Neumann boundary condition, which 
is defined as a special case of Eq. ( |72|) as follows. In 
the c-dimensional space of the vectorial bosonic fields, we 
impose "D" for the c^-dimensional subspace Vd and "N" 
for the remaining djy(= c ~ du)-dimensional subspace V/v 
perpendicular to it. Namely, 

s ■ d x $(x)\T) = for s G V D , (76) 
s ■ d x Q(x)\T) = for s G V N . (77) 

Let P\\ and P± be the projection operators onto Vd and 
Vn, respectively. Then TZ in Eq. ([73) is expressed as 



TZ = IP + (~I)P± = P\ - Pi 



(78) 



which is the reflection operator about the "surface" Vd- 
As explained above, the corresponding boundary state 
|r) is constructed as a linear combination of Ishibashi 
states @: 



gr E I («>#)))> 

(u,v) 



(79) 



where gr is a prcfactor to be determined later and 
the summation runs over all possible (u, v) satisfying 



Eq. (75). Usually, instead of the condition (|75|), it is 
sufficient to require separate conditions for u and v: 



TZu 



-u, IZv 



(80) 



Since u and v live on different lattices S and S, a solution 
(u, v) satisfying only Eq. ( ]75| ) and not Eq. ( p0| ) appears 
only when the primitive vectors of the lattices are fine- 
tuned, and is not considered in the present discussion. 
Because of the definition (|78| ) of 1Z in the present case, 
the conditions (p(I) imply 



u G Vn, v G Vd- 



(81) 



Let Sjv be the set of u G S satisfying u G and 
be the set of v G S satisfying v G Vd- Then, Eq. ( |7^ ) is 
rewritten as 



|r> - 5r E E |(t?,tO)>. 



(82) 



The prcfactorr-ar is fixed by requiring Cardy's consis- 
tency condition.p3 stated as follows. We impose bound- 
ary conditions T and T' at the imaginary time = (3/2 
and respectively, and consider the transition amplitude 
(partition function): 



= <r| 



r') 



(83) 



This can be expressed as a function of 

q = e 27 " T = e- 27r/3/i , (84) 

where r = i/3/L is the modular parameter!! (this picture 
is referred to as the "closed string channel" ). By modular 
transformation, we exchange the roles of space and time 
and express Z^v as a function of q = e~ 2l "/ r = e -27ri /^ 
("open string channel"). In this picture, we may define 
the Hamiltonian Hrr' for a ID system with two bound- 
ary conditions T and V at the ends, and write the parti- 
tion function as Zrr'(q) = Tr e~ LHrl " . This means that 
Zrr 1 (?) is determined by the spectrum of i?rr' • There- 
fore it should have the form 



Z rT ,(q) = J2N£r'Xl h (q), 



(85) 



where Xh"(<l) is a character of the Virasoro algebra. The 
coefficient N r l r , can be interpreted as the number of pri- 
mary fields with conformal weight h, pid has to be a 
non-negative integer (Cardy's conditionE 2 !) . Usually it is 
also required that 7V° r = 1, where h = corresponds to 
the identity operator. This is related to the uniqueness 
of the ground state of i?rr'- This requirement can be 
used to fix <7r- 

Now we calculate the amplitude between two |F)'s de- 
fined by Eq. <M): 



Zrr(q) 



2( 1 
9r 



E E * f(S2+ ^ (86) 
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where 



■ l D" conditions: 



V(q)=q 1/24 lit 1 -?" 1 ) 



(87) 



is the Dedekind 77 function. By modular transformation, 
we can rewrite Zrr using q: 

Z Tr (q) =^ 7 r- c / 2 Wo(H JV )- 1 wo(Hz 3 )- 1 



where vo(...) represents the unit cell volume of the lattice. 
Here we have used the following identities: 



/8 



1/2 



E *** 



1 f /3 A 



.l6 = N 



dw/2 



E ^ 



fl 6Hd 



0_ 



do/2 



(89) 
(90) 

(91) 



The second and third equations come from the multi- 
dimensional generalization of the Poisson summation for- 
mula. To satisfy Cardy's consistency condition above, we 
require the coefficient of the term with (r, s) = (0, 0) to 
be unity, obtaining 



3r 



(92) 



The constant g^ appears in the overlap between the 
ground state jvP) of if and the boundary state (82): 
gr = (^|r). This means that the partition function Zrr 1 
has a multiplicative constant contribution <?r5r' com- 
ing from the boundaries in the limit j3/2 3> L ^> 1. 
This result can be interpreted as follows. In the open 
string channel picture, the ID system described by Hrr' 
has the "spacial length" (3/2 and the "inverse tempera- 
ture" L. The ground state of i?rr' is unique, and there- 
fore the thermal entropy goes to zero in the "zero tem- 
perature" limit 1/L — s- 0. On the other hand, when 
/3/2 ^> L ^> 1, the "temperature" 1/L is high enough 
and the spectrum of Htv looks effectively continuous. In 
this case, the thermal entropy acquires a constant contri- 
bution log(gr.<?r')j i n addition to the standard extensive 
contribution linear in temperature. Because of this, gr is 
referred to as the boundary ,-iiground-state degeneracy," 
and is generally non-integer .cil 



D. Boundary conditions Ti and T2 

The two boundary conditions I\2 in Eqs. ( |59| ) and 
( |60| ) can be expressed as special cases of "mixed" Dirich- 
let /Neumann conditions discussed in the previous subsec- 
tion. The condition T2 in Eq. (|6C]) leads to the following 



- -ft 

r 

'(5) = TXV\ 



Q = d x {^-' l -^>) = -_{b 2 



6? j) )-&$ 



(93) 



= 1,2; br=b?>). 

Therefore, the subspace Vd2 where "D" is imposed is 
spanned by the following (non-orthogonal) basis vectors: 



s x = b y 2 > -b y 2 >, s 2 = b\*>-b 



m _ r(3) 
1 u i > 



S3 — 2 — 2 , S4 — o 1 — o 1 



(94) 



In the perpendicular space Vn 2 , the field $ is free at the 
boundary, and therefore we impose "N" , where the dual 
field O is locked. This space is spanned by 



ta 



7<D 



?(3) 



? (2) 
*2 5 



?(2) 



ti = a\ 



(4) 



.(3) 
l l > 

f (l) 



(95) 



Now we consider the lattices Sjv2 and ^D2 used to con- 
struct the boundary state |r 2 ) as in Eq. (p2]). Since any 
linear combination of {sj} with integer coefficients be- 
longs to the lattice 3, {sj} can be used as the primitive 
vectors of S£> 2 (= 3 fl Vd2)- Similarly, {t,} can be used 
as the primitive vectors of 5^2 (= 3 fl Vn2)- 

What are the meanings of the lattices 5jv2 and 3^2 
introduced in this way? Initially, in the mode expansions 
(|63|), the eigenvalue u of u can take any element of the 
lattice 3 and thus be expressed as 



EE 4m 

J = l v = \ 



(96) 



(i) 

where n), is an integer representing the winding number 
of 4>v in Fig. ^(a) in the x direction. After imposing 
the boundary condition T2, u lives on the reduced lat- 
tice S7V2, where as indicated by Eq. (|95|), the winding 
numbers obey the constraints: 



(1) (2) 



(2) (3) 



(3) 



,( 4 ) 



(97) 

"2 — "2 3 n i ' — n i ' ■ 
These equations simply mean that the winding numbers 
of (/)„ 's on two sheets connected through T2 in Fig. |^(a) 



should take the same integer. Similarly, Eq. (|94|) implies 

that the winding numbers m^' of 9v's on two sheets 
connected through T2 should take mutually opposite in- 
tegers. 

Using the primitive vectors {s*/}, the unit cell volume 
of 3c2 is calculated as 



51 ■ Si S\ ■ s 2 Si ■ S3 Si • s 4 

5 2 • S*l S 2 • S 2 S 2 • S3 S 2 • S4 

53 • Si S3 • S 2 S3 • S3 S3 • S 4 

54 • Si S4 • S*2 S4 • S3 S4 • S4 

det ( 2f 2 K] /2 K 1 J 2 M 2 



(98) 
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where M2 is the 4x4 matrix defined in Eq. (p7|). Since 
the general case of integer n > 1 can be handled by sim- 
ply replacing M 2 by the 2n x 2n matrix M n [defined in 
Eq. (p7|)1 , we proceed our discussion in this general case. 
Now the boundary conditions I\2 are imposed on a 4n- 
component boson. We obtain 



v (~ D2 )=(2^K 1 J 2 K 1 J 2 



(detM„) 1/2 . (99) 
Similarly, we obtain the unit cell volume of Sjv2 as 

v (E N2 ) = (2r 2 K- 1/2 KZ 1/2 ) n (detM„) x / 2 . (100) 
Therefore, using Eq. (|92|), the factor <7p 2 i s calculated as 
gr 2 = ^ l i>o(S D2 ) 1 / 2 z; (S J v 2 ) 1 / 2 = (det M^ 1 / 2 . (101) 
A similar procedure for T\ yields <?n = 1. 

E. Calculation of reduced density matrix moments 

We consider the transition amplitudes, Zt 1 v 2 an d 
Zy^x- The calculation of Zr 1 r 2 f° r arbitrary f3 is a dif- 
ficult issue because the 1Z matrices for the two boundary 
conditions do not commute with each other. However, 



as mentioned at the end of Sec. IV C , one can still derive 
the asymptotic expressions in the limit ^>I>1 (i.e., 
q — > 0). The results are 



Z Tl T 2 W (r i |*)g- 4 ™/24 (x j / | r2) 
Zr lFl « (r i |*)g- 4 «/24 (x j / | ri) 
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-4n/24 



<? ri <?r 2 , (102) 



g" 4 "/ 24 ^,, (103) 



from which we obtain 



Tip* 



ffr 2 

0Ti 



(detM„) 1 / 2 . 



(104) 



This constant is exactly the same with that appearing 
in Eq. (|3l]). So far, we have been concerned only with 
the regulated part of Tr p\ and have neglected diver- 
gent contributions from the short-range cutoff. In gen- 
eral, the logarithm of the partition function, log Zrv , 
for a cylinder contains terms proportional to the area 
t;L and the circumference L (Refs. |2^j55| ). The coef- 
ficient of the circumference term depends on the de- 
tails of the boundary conditions while that of the area 
term depends only on the bulk properties. Therefore, in 
— log(Tr Pa) = — log(Zr 1 r 2 /Zr 1 r 1 ), the area terms can- 
cel out while the circumference terms do not, leaving a 
contribution aL. In this way, the lin ear co ntribution in 
S n (with integer n > 2) found in Sec. [II C is also repro- 
duced. 



NUMERICAL ANALYSIS 



In this section, we test the analytical predictions of 



0.5 




0.1 



l/L 



0.2 



FIG. 



5: (Color online) 1/K±(L) [Eq. ( |l08| )] versus l/L for 
po = 1/3 and V = — 1. Filled and empty symbols show the 
data of 1/K+(L) and l/K-{L) (with L = 6,9,12,15), re- 
spectively. Lines show the fitting with the quadratic form 
1/K±{L) = l/K± + a/L + b/L 2 . Our motivation to plot 
1/K±(L) instead of K±(L) is that the formers vary in a 
smaller range [0, 2] in the parameter range of our interest. 



a hard-core bosonic model on a ladder. The Hamiltonian 
of the ladder model is given by 



ff=EE -t(b] >u b j+1 , v +h. c ) 

„=l,2j=l L 



(105) 



Sees. [II and IV in a numerical diagonalization analysis of 



where bj tl/ is a bosonic annihilation operator at the site 

j on the v-th leg, and rij tV = bj v bj, v is the number 
operator defined from it. Here, t and V represent the 
hopping amplitude and the interaction between nearest- 
neighbor sites on each leg, and U represents the interac- 
tion along a rung. We impose the hard-core constraint 
tij v = (bj „) 2 = 0, and therefore the bosonic operators are 
equivalent to spin-i operators as bj V = S~^, bj v = S^ v . 
We assume the PBC &z+i,t/ = b\ iV . We define the aver- 
age particle density as p = (Ni + N 2 )/(2L), where N v 
is the particle number on the ^-th leg. We assume t > 0, 
—2 < V < 0, and U > 0; this case was studied recently 
in Ref. We set t = 1 in the following. As explained in 
Appendix [b| and in Ref. [56[ this model is equivalent to 
a fcrmionic model on a ladder under the Jordan- Wigner 
transformation. In particular, for V = 0, the model is 
equivalent to the S'J7(2)-symmetric fermionic Hubbard 
chain, which is solvable by Bethe ansatz. In the Hubbard 
chain, the two legs u = 1 , 2 are identified with the spin- 
up/down states, and the symmetric/antisymmetric sec- 
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1 /K + • 

1 /K_ o 




(a) 



FIG. 6: 1/K± for po — 1/3 and V = —1, obtained by extrapo- 
lating finite-size data as in Fig. ||. Lines show the perturbative 
estimates (Jic|) . 



tors correspond to charge and spin modes, respectively. 

We briefl y review the recent results of Rcf. |5(] on the 
model (105). For U — 0, the model decouples into two in- 
dependent Bose gases, each equivalent to a solvable spin- 
i XXZ chain in a magnetic field. Each chain forms a TLL 
described by the Hamiltonian (|^). The velocity v and the 
TLL parameter K of . ea ch XXZ chain can be determined 
from Bethe ansatz.EHH For small U > and po 1/2, 
the inter-chain coupling can be analyzed along the same 
argument as Sec. II A, leading to the perturbative esti- 
mates ( |io| ) of the rcnormalized velocities v± and TLL pa- 
rameters K± (here, the lattice constant is set to unity). 
As seen in this estimate, K_ increases with increasing 
U . For V < 0, it was found that K- finally diverges 
as U approaches certain U c , where a first-order phase 
transition to a population-imbalanced state (Ni ^ N2) 
occurs. Here we focus on the uniform phase {N\ = N2) 
in < U < U c described by the effective Hamiltonian 
in Eqs. (||) and In the solvable case V = 0, the 

transition is known not to occur, and the uniform phase 
continues for arbitrary large U > 0. In our calculation 
presented below, we fixed the density at po = 1/3, and 
performed calculations for V = —1, —0.5, and 0. 



Before presenting our results on entanglement, let us 
explain our method for calculating the TLL parame- 
ters K±. In the solvable case V = 0, K± can be de- 
termined accurately by numerically solviag-thc integral 
equations obtained from Bethe ansatz.&O For other 
cases, we determined K± in numerical diagonalization 
of finite systems (up to L = 15) by using the method of 
Rcfs. 62 63|. In this method, we define fij t ± := fij t \ in2j 
y — po, and examine their correlation func- 
[hj t ±hj+r t ±) . Using the bosonic represen- 



with hjjy := r 
tions C± (r) := 

tation of operators, these correlation functions arc shown 
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FIG. 7: (Color online) Si and Soo versus L for po = 1/3 and 
V = — 1. Solid lines show the fits with the scaling form (10£) 
and broken lines show the linear part a„L + 7„. 



to have the asymptotic forms 
K ± A 



C±(r) 



± 



(7rr) 5 



cos(2kpr) + 



(106) 



where kp := 7rpo is the Fermi momentum in the cor- 
responding fcrmionic model and A± are non-universal 
coefficients. In the SU (2)-symmctric case, a marginally 
irrelevant perturbation produces rnxdimlicative logarith- 
mic corrections in the second term.EltJO Performing the 
Fourier transform, only the first term contribute for a 
small wave vector q, leading to 



-iqr 



-\q\ (}»0). (107) 



In a periodic finite-size system of length L, we evaluate 
this for q = 2n/L, leading to the finite-size estimate of 
the TLL parameters: 



(108) 
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FIG. 8: (Color online) Si versus L for p = 1/3 and V = 0. 
Solid lines show the fits with the linear form S„ = a n L + 7„. 



The data of 1/K±(L) are extrapolated into L —> oo as 
illustrated in Fig. g. The values of l/if± obtained by 
the extrapolation are plotted as a function of U in Fig. ||, 
in reasonable agreement with the perturbative estimates 
@ for small U(< 1). 

Let us now present our results on the (Renyi) entan- 
glement entropies S n (with n = 1, 2, oo) between the two 
legs of the ladder. These entropies are calculated in the 
ground states of finite-size systems (up to L = 12) ob- 
tained by Lanczos diagonalization. The data of S n well 
obey a linear function qL,L. For V = — 1 and —0.5, we 
find that a scaling formL3 



5/i 



(109) 



fits the data very well as shown in Fig. |7| The linear 
part a n L + j n (broken lines) crosses zero around L = 3, 
whic h mea ns that the short-range cutoff eto discussed in 
Sec. IIIB is given by eto ~ 3. For V = 0, in contrast, 
a simple linear form S n = a n L + 7„ fits the data better 
as shown in Fig. || The extracted constant 7„ is plot- 
ted as a function of U in Fig. [| Using the values of 
the TLL parameters K± obtained numerically, the for- 
mulae of 72 and 7oo in Eqs. (|34|) and ( |35| ) are also plotted. 
For V = — 1 and —0.5 [FigTg(a), (b)], we find a broad 
agreement between the numerical data and the analytical 
formulae. The difference between them are within sa 30% 
of their values. We note that our calculations of both j n 
and K± are based on finite-size systems with L < 15. 
We expect that calculations in larger systems (by using, 
e.g., the quantum Monte Carlo method of Ref. g6|) would 
demonstrate a more accurate agreement with the analyt- 
ical predictions. For V = (the Hubbard chain case), on 
the other hand, we find a significant difference between 
the numerical and analytical results — the numerical re- 
sults are roughly four times as large as the analytical 
results. The origin of this significant difference occurring 
only for V — will be discussed later in this section. 
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FIG. 9: (Color online) 7„ (with n — 1,2, cxd) as a function of 
U , for po = 1/3 and V = —1, —0.5, 0. The analytical formulae 
of 72 and 700 in Eqs. (^) and ( |3q ) are also plotted using the 
values of the TLL parameters K± obtained numerically. 



In Fig. |To|(a) , we plot the relation of 7„ and K+ / K- 
using the data for V = — 1 and V — —0.5. We can 
again confirm that for 72 and 700, the numerical data 
and the analytical formulae show a broad agreement. 
Furthermore, we observe that the data of 71 for two val- 
ues of V show a broad agreement, which suggests a uni- 
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One can see in Fig. 10(b) that the data for n — 1.2 and 
1.4 indeed intervene between the data of n = 1 and 2. 
The issue of how the small- re behavior of 7„ changes in 
the range 1 < n < 2 is subtle within the present data. We 
here propose two possible scenarios and leave the issue 
open for future studies. One scenario is that the exponent 
b decreases smoothly in the range 1 < n < 2 although it 
is fixed at b — 2 for n > 2. Another scenario is that the 
quadratic behavior of Eq. fl45[ ) holds for arbitrary n > 1, 
but the range of re where the quadratic term dominates 
shrinks gradually as n approaches 1. 

Finally, let us discuss the origin of the significant differ- 
ence between numerical and analytical results observed 
for V = (the Hubbard chain case) in Fig. |^(c). In 
5L/(2)-symmctric systems like the Hubbard chain, it is 
known that a marginally irrelevant perturbation pro- 
duces non-trivial corrections to the predictions of the 
pure Gaussian model in various physical quantities. In 
particular, its effects are enhanced in the presence of non- 
trivial boundary j^onditions, as discussed in thn_spin-i 
Heisenberg chainElrGJ and the Hubbard chainEjO with 
open ends. In the present case, the system has a sim- 
ple periodic boundary condition in space, and non-trivial 
boundary conditions are imposed in the imaginary time 
direction as presented in Sec. [V. We expect that a 



0.05 



0.1 



0.5 



pcrturbativc calculation using boundary states, as was 
done in Rcf. |6^| would clarify non-trivial effects of the 
marginally irrelevant perturbation. 



FIG. 10: (Color online) (a) 7„ versus K+/K-. (b) -7„/re 2 
versus re := (K- — K+)/(K- + K+) in logarithmic scales. 
The data of 71 are fitted with the form (^). The analytical 
formulae of 72 and 700 are from Eqs. (B4) and (J35j) . 



versal relation between 71 and K+/K_ 
we have expanded 7„ (with n — 2, 3, . 



:= (K_ -K + )/{K_ 



In Sec. [Ill q , 

) in terms of 
and found the leading 



dependence ( 

,.2 



Motivated by this observation, we 
plot -~7„/k 2 as a function of re in logarithmic scales 
in Fig. |o|(b). This figure also presents some data for 
1 < n < 2 obtained in a similar way. As expected, the 
data for n — 2 and 00 stay around constants as re de- 
creases, although these constants arc slightly larger than 
those expected from Eq. ([l5]). The data for n — 1, how- 
ever, increase as ft decreases, and follow straight lines 
in logarithmic scales in Fig |l^(b) . We fit the data with 
the form 71 = —are 6 (as mentioned in Eq. (|l6|)) in the 
range < ft < 0.5, obtaining (a, b) w (1.13,1.70) and 
(1.16,1.62) for V = — 1 and 0.5, respectively. This in- 
dicates that the leading ^-dependence of 71 contains a 
non-trivial exponent b ~ 1.6-1.7, in marked contrast to 
the quadratic dependence (0) of 7„ for integer n > 2. 

In spite of the qualitatively different small-re behaviors 
for n = 1 and n > 2, we have found that for fixed re, 7„ 
changes rather smoothly when n is changed from 2 to 1. 



VI. SUMMARY AND DISCUSSIONS 

We have considered two coupled TLLs on parallel 
chains and calculated the Rcnyi entanglement entropy 
S n between the two chains. Wc formulated the problem 
in the path integral formalism, and related S n with inte- 
ger n > 2 to the partition functions on certain non-trivial 
manifolds. These partition functions were calculated us- 
ing two analytical methods. We argued that S n obeys a 
linear function of the chain length L followed by a uni- 
versal sublcading constant j n . The two methods led to 
the same formulae for j n , which are written as functions 
of the ratio of TLL parameters. The obtained formulae 
were checked numerically in a hard-core bosonic model on 
a ladder. When the model is away from the SU(2) case, 
the numerical data of 72 and 700 showed a broad agree- 
ment with analytical formulae. The agreement among 
two analytical approaches and numerical results has of- 
fered a convincing evidence of the universality of j n with 
integer n > 2. Our numerical results also suggested that 
the sublcading constant 71 in Si is also universal and that 
its leading dependence on re := (K- — K + )/{K^ + K + ) 
obeys a non-trivial power function, in contrast to the 
quadratic dependence of j n for integer n > 2. In the 
SU (2)-symmetric case, the numerical data of 72 and j n 
differ significantly from the analytical formulae, which 
indicates a strong effect of a marginally irrelevant per- 
turbation. 
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Recently, it has been discussed that the particle num- 
ber fluctuations in a subsystem show similar scaling 
behavior— to the entanglement entropy in a number of 
systemsJlj This is an interesting proposal relating the 
entanglement entropy to an experimentally observable 
quantity. In our setting of two coupled TLLs, parti- 
cle number fluctuations in a chain are completely absent 
since the particle number is separately conserved in each 
chain. On the other hand, finite entanglement entropy 
does exist between the two chains, and obeys a linear 
scaling with the chain length L as we have discussed. 
Therefore, our study offers a counterexample to the sim- 
ilarity of the two quantities. We comment that different 
behaviors of the two quantities have also been discussed 
in the dynamics of fractional quantum Hall states after 
a local quantum quench.E3 

Our formulations for studying two coupled TLLs can 
be extended to study the entanglement in multicompo- 
nent TLLs. An exciting possibility is to study the. entan- 
glement entropy in a sliding Luttinger liquidjxu which 
appears in a 2D array of coupled TLLs. To be specific, 
we define such a system on a torus of length L x and L y 
in two directions. Here, TLLs, described by the bosonic 
fields 4>j{x) with j = 1, 2, . . . , L y , are running along the 
x direction and are mutually coupled in the y direction. 
Assuming the translational invariance in the y direction, 
it is natural to introduce the Fourier transform of the 
bosonic fields in the y directions: 



x 2 - x\ 



L x /2, we predict a scaling 



4> q (x) = 



(110) 



log 



7T 



ir(x 2 - x-i) 



const. 



(Ill) 



= — \ogL x + const. 



In these ways, the entanglement entropy shows qualita- 
tively different scaling behaviors depending on in which 
direction one cuts the system. Such a highly anisotropic 
character of entanglement is related to the anisotropic 
correlations in this systcpi,|-ajid is in marked contrast to 
non-interacting fcrmionalxllj and Fermi liquids. 113 
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Appendix A: Ground state wave functional of a TLL 



with q = 2im y /L v [n y = 0, 1, . . . , L y — 1). In a sliding 
Luttinger liquid, the total Hamiltonian decouples into 
independent TLLs, each defined for <f> q with the renor- 
malized TLL parameter K q and the velocity v q . Now 
we consider dividing the torus into two cylinders of the 
same size by cutting it along two lines either in the x 
or y direction. Cutting along x is similar to the prob- 
lem of this paper; it can be treated by generalizing the 
formulation in Sec. |y| using more complicated "mixed" 
Dirichlet/Neumann boundary conditions. It then leads 
to the linear scaling of the entanglement entropy with L x , 
followed by a subleading constant determined by L y TLL 
parameters. The coefficient of the linear term can depend 
on L y , but we expect that it converges to a constant for 
sufficiently large L y because of the short-range character 
of the correlations in the y directions. When we cut the 
system along y, the original bosonic fields <j>ji x ) are cu t 
at the same positions (say, x = x\ and x 2 ) independent 
of j. Then the Fourier components (j) q (x) are also cut 
at the same positions for all <?'s. Therefore, the entan- 
glement entropy in this case can be treated in the same 
way as the single-interval entanglement entropy in a ID 
gaplcss system with central charge c =r L~. Using the 
finite-system formula in the latter caseMli3 and setting 



Here we consider a single-component free boson Hamil- 
tonian (defined by the fields <p(x) and 9(x)) with the 
TLL parameter K, and derive the expression of the 
ground-state wave functional (<p\^)- Such wave fuj, 

path integral 




and the Calogero- 
This problem is also 



als have been derived by us: 
the Schrodinger formalism 
Sutherland wave function.! 
closely related to the effective action for the boundary 
degrees otfrcedom discussed in the cantext of dissipation 
problemsLa and impurity problems E3 Here we present a 
simple derivation in the operator formalism. Since the 
winding numbers (zero modes) of the bosonic fields are 
zero in the ground state, we ignore them in the following 
discussion. 

The field <fi is expanded as 



m = E 



47rm 



with k m = 2Ttm/L and 



[ tt m, a t'l 



ik m x~\ 



(Al) 



= a 



r ' T l - f) 
i'J — °i> 



This is a one-component version of Eq. (p3(). The ground 
state |*) is defined by a m /R \^) = (Vm g N). By anal- 
ogy with the quantum mechanics of a harmonic oscillator, 
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we introduce the "coordinate" and "momentum" opera- Then the normalization factor j\f is calculated as 



tors, Xm R and Pm R , for each mode via 



X L I R 4- 4 P L / R 

71 ' 



yL/R _ -pL/R 

7T~ 



(A2) 

The hermittian operators Xm R and Pm R satisfy 
the canonical commutation relations [X^P^,] — 
[X R ,P R ,] = iS mm f. We further introduce 



pR _j_ pL 



yn _i_ y i. 

m,± = f= , P m ,± = 7= j ("1 > 0) 



V2 



V2 



(A3) 

which are related to the "center of mass" and "relative" 
motions of the left /right-moving modes labeled by m. 
Then, Eq. ( |Al| ) is rewritten as 



o = E 



K 



(X ro ,+ +iP m ,-)e l 



+ (X m>+ - iP m _)e 



-ik m x 



(A4) 



This expression "diagonalizes" </>(x) because all A mj+ 's 
and P m> _'s commute with each other. 
The state |y>) is defined by 



(x)\<p) = <p{x)\tp) (0<x<L). 



(A5) 



Prom Eq. (A4), one can see that \tp) is given by a simul- 
taneous eigenstate of {X m> y, P TOi _} to> q. We expand the 
field configuration y(a;) as 



9 



1 OO 

VJj m=l 



~ Ak m x i ~* „ — ik m x\ 

ip m e +f m e m J 



(A6) 



Then the coefficient <^ m is related to the eigenvalues, 
A m-+ and P m _, of X m + and P m _ as 



A m:+ + 1P TO _ = -y/ —j^$r, 



(A7) 



From the solution of a harmonic oscillator, the ground 
state wave function is written in a Gaussian form in terms 
of Xm + 's and P m _'s as 



({A mi+ ;P m ,_}|*) oc cxp 



OO 



^2 , p2 

m— 1 



(A8) 



The wave function in terms of (p m 's is then given by 

1 / 1 °° \ 

(#U|*) = ^cxp I -— J2 kmWm?\ ■ (A9) 

We normalize the wave function such that 

/ H(dv m d<p* m ) |({^ m }|*)| 2 = l. (A10) 



OO p 

Af = ] / d{p m dtp* m exp 

m=l ^ 



2k 



K 



mi. |2 



It is interesting to transform Eq. (| 
space representation: 



(vl*> 



1 



m— 1 

(All) 
into the real- 



(A12) 



with 

S[<p] 



L pL 

dx\ I dx2 



1 

2^ Jo 

d x ip(x 1 )d x (p(x 2 )loi 



e L 



(A13) 



Using the charge density measured from the average, 
Sp(x) = —d x ip(x)/y^TT, this is rewritten as 



1 



8p(x\)dxi I 5p(x2)dx2 log 



e e 



(A14) 

This can be viewed as the energy of a classical Coulomb 
gas placed on a unit circle with a logarithmic repulsive 
potential. Such a Coulomb gas structure of the ground 
state wave function is directly seen in the Jastm\K=jtypc 
ground states of the Calogeipr&utherland modeltErlEl and 
the Haldane-Shastry modelEaS. More detailed discus- 
sions on these connections can be found in Refs. ^3 36. 



Appendix B: Jordan- Wigner transformation for a 
ladder 



Under the Jordan- Wign er tr ansformation, the hard- 
core bosonic model in Eq. (105) is equivalent to a spin- 
less fermionic model on a ladder, where all the bosonic 
operators b^ v in Eq. ( p_05| ) are replaced by fermionic ones 
fj jU . This transformation is defined as 



fj,i = ex P 
fj,2 = exp 



i=i 

/ L J-l > 

in y^Ttj.i +^Jii,2 
\i=i 1=1 / 



(Bl) 

bj,2, (B2) 



where the "string" part runs first along the first leg and 
then al ong the second leg. In particular, for V = 0, the 
model (105) is equivalent to the solvable fermionic Hub- 
bard chain, where the two legs v = 1, 2 are identified 
with the spin- up/down states. Although the Hamilto- 
nian retains the same form under this transformation, the 
boundary condition is transformed in a non-trivial way. 
For example, the PBC bh+i,v = &i i/ on the bosons cor- 
responds to the boundary condition /l+i,v = & % fi,v 
on the fermions, where N v is the number of particles on 
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the u-t h leg . Our motivation to consider the bosonic 
model (105) instead of the fermionic one is that in the 
unifo rm phase which we consider here, the bosonic model 
(105) with the PBC has a unique ground state, irrespec- 
tive of the chain length L and the total particle num- 



ber N = Ni + N 2 . On the other hand, for {7 = 0, the 
fermionic model with the PBC has degenerate ground 
states for some L and N. Although this degeneracy is 
split for U > 0, some irregular size dependence occurs as 
a remnant of the degeneracy at U = 0. 
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